(12) INTERNATIONAL APPLICATION PUBLISHED UNDER THE PATENT COOPERATION TREATY (PCT) 



(19) World Intellectual Property 
Organization 

International Bureau 

(43) International Publication Date 
7 October 2004 (07.10.2004) 




PCT 





II 



(10) International Publication Number 

WO 2004/086310 Al 



(51) International Patent Classification 7 : 



G06T 17/00 



(21) International Application Number: 

PCT/AU2004/000342 

(22) International Filing Date: 19 March 2004 (19.03.2004) 

(25) Filing Language: English 

(26) Publication Language: English 



(30) Priority Data: 

2003901625 



28 March 2003 (28.03.2003) AU 



(71) Applicant (for all designated States except US): COM- 
MONWEALTH SCIENTIFIC AND INDUSTRIAL 
RESEARCH ORGANISATION [AU/AU]; Limestone 
Avenue, Campbell, ACT 2612 (AU). 

(72) Inventors; and 

(75) Inventors/Applicants (for US only): LI, Rongxin 
[CN/AU]; U/129B Park Road, Dundas, NSW 2117 (AU). 
OURSELIN, Sebastien [FR/AUJ; 30 Parkview Road, 
Fairlight, NSW 2094 (AU). 

(74) Agent: SHELSTON IP; 60 Margaret Street, Sydney NSW 
2000 (AU). 

(81) Designated States (unless otherwise indicated, for every 
kind of national protection available): AE, AG, AL, AM, 



AT, AU, AZ, BA, BB, BG, BR, BW, BY. BZ, CA, CM, CN, 
CO, CR, CU, CZ, DE, OK, DM, DZ, EC, EE, EG, ES, FI, 
GB, GD, GE, GH, GM, HR, HU, ID, IL, IN, IS, JP, KE, 
KG, KP, KR, KZ, LC, LK, LR, LS, LT, LU, LV, MA, MD, 
MG, MK, MN, MW, MX, MZ, NA, NI, NO, NZ, OM, PG, 
PH, PL, PT, RO, RU, SC, SD, SE, SG, SK, SL, S Y, TJ, TM, 
TN, TR, TT, TZ, UA, UG, US, UZ, VC, VN, YU, ZA, ZM, 
ZW. 

(84) Designated States (unless otherwise indicate^, for every 
kind of regional protection available): ARIPO (BW, GH, 
GM, KE, LS, MW, MZ, SD, SL, SZ, TZ, UG, ZM, ZW), 
Eurasian (AM, AZ, BY, KG, KZ, MD, RU, TJ, TM), Euro- 
pean (AT, BE, BG, CH, CY, CZ, DE, DK, EE, ES, FI, FR, 
GB, GR, HU, IE, IT, LU, MC, NL, PL, PT, RO, SE, SI, SK, 
TR), OAPI (BF, BJ, CF, CG, CI, CM, GA, GN, GQ, GW, 
ML, MR, NE, SN, TD, TG). 

Declaration under Rule 4.17: 

— of inventorship ( Rule 4. 1 7( iv)) for US only 

Published: 

— with international search report 

For two-letter codes and other abbreviations, refer to the "Guid- 
ance Notes on Codes and Abbrexdaiions" appearing at the begin- 
ning of each regular issue of the PCT Gazette. 



(54) Title: COMBINING FRONT PROPAGATION WITH SHAPE KNOWLEDGE FOR ACCURATE CURVILINEAR MODEL- 
LING 



10 



< 



00 




11 



12 



Construct Propagation Map 



13 



(57) Abstract: A method of constructing a curvilinear model 
of a tubular structure from a 3dimensional data set, the method 
comprising the steps of: (a) forming a substantially minimal cost 
path for the tubular structure by the steps of: (al) constructing 
a distance field; (a2) utilising the distance field to construct a 
propagation map; (a3) utilising the propagation map to construct 
the substantially minimal cost path; (b) for each point on the sub- 
stantially minimal cost path, providing the steps of: (bl) sampling 
the 3dimensional data set in a substantially locally perpendicular 
axis to the minimal cost path; (b2) constructing a deformable tube 
around the substantially minimal cost path determined by alter- 
ations in pixel values around the minimal cost path; (b3) extract- 
ing an axial measure of the deformable tube as the curvilinear 
model of the tubular structure. The method can further include 
the step (M) of transforming the curvilinear model back into the 
original data set. 
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COMBINING FRONT PROPAGATION WITH SHAPE KNOWLEDGE FOR 

ACCURATE CURVILINEAR MODELLING 

FIELD OF THE INVENTION 
The present invention relates to the field of automated extraction processes in 
5 computer graphics or the like. In particular, the present invention discloses the 

automated extraction of curvilinear models of tubular structures in 3D data sets such as a 
series of scanned medical images or the like. 

BACKGROUND OF THE INVENTION 
In neurosurgery or endovascular surgery, vessel system understanding of the 
10 patient is of fundamental importance. Further, accurate curvilinear modelling of 

anatomical tubular structures may become a crucial step in computer- assisted surgery 
planning. Curvilinear models of a tubular entity is an abstraction in object 
representation that retains only those geometric properties that are essential to length, 
angle and tortuosity measurements while excluding all potentially interfering details 

» 

15 such as widths, width changes (e.g. aneurysms) and branching. Previous methods for 
extracting centrelines (or medial axes) have in many situations not been able to construct 
an accurate curvilinear model. 

A classical method for curvilinear modelling is the skeletonisation technique. This 
is often achieved by topology-preserving thinning using an algorithm that deletes simple 

20 cells after region growing ( see N. Nikolaidis and I. Pitas. 3-D Image Processing 
Algorithms. John Wiley, 2001). 

Alternatively voxels in the middle of a curvilinear structure can be located by 
smoothing the image using a multi-scale response. A model-based scale-space filtering 
approach has been also proposed by several authors. The centreline is located by 
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detecting optimal responses derived from Hessian matrices at a range of scales of 
interest. This class of algorithms show good results but some local problems occur at 
junctions or tangent vessels. 

Direct tracking from a seed point has also been used to extract the optimal path. 
5 From a given point, the most likely local orientation of the vessel can be searched or 
predicted. This involves an iterative search, prediction, centering, and image re- 
sampling. However, the fidelity of the model to the "true" centerline thus extracted 
depends heavily on the centering algorithm. 

Graph-search principles have also been used to calculate the optimal path. To 
10 achieve this task, Dijkstra's algorithm has been the most efficient and widely used one 
for this purpose. 

A recent approach proposed by T. Deschamps and L.D. Cohen "Fast extraction of 
minimal path in 3D images and applications to virtual endoscopy", Medical Image 
Analysis, 5: 281-299, 2001 incorporates multiple improvements over the standard fast 

15 marching algorithm, including a technique of using multiple passes of fast marching 
(from different initial positions) to centre the track. 

Unfortunately, using existing centreline extraction algorithms, a high level of 
accuracy cannot be guaranteed. A common drawback of the existing techniques is they 
are not sufficiently resistant to neighbourhood interference. Examples of the types of 

20 neighbouring interference are shown schematically in Fig. 1 with an initial Vessel 1 
having a centreline 2. The vessel and centreline is then shown distorted by a branch 4, 
another anatomical structure 5, vessel compression 6 and by an aneurysm 7. 
Consequently, existing algorithms may introduce spurious tortuosity at places of 
asymmetricity, e.g. vessel branching, aneurysms, neighbouring objects touching the 
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object of interest. This may significantly impact the accuracy of any length and 
tortuosity measurements based on the model. 

SUMMARY OF THE INVENTION 
It is an object of the present invention to provide for an improved method of 
5 curvilinear modelling of tubular type structures. 

In accordance with a first aspect of the present invention, there is provided a 
method of constructing a curvilinear model of a tubular structure from a 3-dimensional 
data set, the method comprising the steps of: (a) forming a substantially minimal cost 
path for the tubular structure by the steps of: (al) constructing a distance field; (a2) 
10 utilising the distance field to construct a propagation map; (a3) utilising the propagation 
map to construct the substantially minimal cost path; (b) for each point on the 
substantially minimal cost path, providing the steps of: (bl) sampling the 3 -dimensional 
data set in a substantially locally perpendicular axis to the minimal cost path; (b2) 
constructing a deformable tube around the substantially minimal cost path determined by 
15 alterations in pixel values around the minimal cost path; (b3) extracting an axial measure 
of the deformable tube as the curvilinear model of the tubular structure. The method can 
further include the step (b4) of transforming the curvilinear model back into the original 
data set. 

The step (b2) further can comprise the steps of constructing a resample data 
20 around the substantially minimal cost path and then constructing the deformable tube 
within the resample data. The distance field can be constructed utilising a Chamfer 
algorithm for a distance measure. The distance field can be thresholded. The propagation 
map can be constructed utilising Fast Marching Method techniques. 
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The tubular structure can comprise a tubular structure within the human body and 
the 3-dimensional data set can comprise a scan of a portion of the human body. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Preferred embodiments of the present invention will now be described with reference to 
5 the accompanying drawings in which: 

Fig. 1 illustrates schematically various forms of possible interference present inside the 
human body; 

Fig. 2 is a flow chart of the process of extraction of a minimum cost path of the preferred 
embodiment; 

10 Fig. 3 is a flow chart of the steps of formation of a deformable tube of the preferred 
embodiment; and 

Fig. 4 illustrates one form of computer system suitable for implementing the preferred 
embodiment. 

DESCRIPTION OF PREFERRED AND OTHER EMBODIMENTS 
15 The method of the preferred embodiment for accurate curvilinear modelling 

substantially reduces the problems described earlier with the existing approaches by 
using shape knowledge in addition to wave propagation. This is provided by a modified 
tubular deformable model to implement the constraints. 

The preferred, embodiment is described with reference to the person skilled in the 
20 art of interactive computer graphics design and medical image processing. In particular, 
it is assumed that the person skilled in the art will be readily familiar with the techniques 
discussed in papers such as T. Deschamps and L.D. Cohen, "Fast extraction of minimal 
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path in 3D images and applications to virtual endoscopy" Medical Image Analysis , 5: 
281-299, 2001 . 

Fig. 2 illustrates a flow chart in the initial steps of the preferred embodiment in 
extraction of a suitable minimal cost path. A distance field is first constructed 1 1 and 
5 then uniformly thresholded 12 to obtain suitable propagation channels. A propagation 
map is then constructed 13 from a start point using Fast Marching Methods (FMM), and 
then backtracked when the end is reached to generate the minimal cost path (MCP) 14. 

Fig. 3 illustrates the subsequent processing steps 20. At each point on the MCP, 
the original data is resampled 21 in planes perpendicular to the local orientation of the 

10 MCP. A deformable tubular model 22 is constructed inside the resampled tube and 
allowed to deform in order to recover the object of interest. Finally, the curvilinear 
model is obtained by extracting the spine of the tubular model and transforming it back 
into the original data. This process ensures symmetric rapid changes, such as a sharp 
turns, are not lessened while many asymmetric changes, such as a branching event or an 

15 aneurysm, can be removed to a large extent. In addition, deforming the model in the 
transformed data domain vastly simplifies the computation, making it efficient. 

Instead of being used to extract the centreline through topology-preserving 
thinning or ridge enhancement as provided previously, distance transforms are used in 
the method of the preferred embodiment as a pre-processing step for larger structures to 

20 reduce the widths of propagation channels. Doing so achieves two purposes. Firstly, 
this helps the tubular model in the next stage to maintain faithfulness to the structure of 
interest, especially at high curvature parts of the structure, and promotes faster 
convergence of the deformation of the model. This is because better initialization for the 
tubular model can be achieved with narrower propagation channels relative to the axial 

25 curvature. Secondly, higher efficiency can be achieved by reducing the size of 
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propagation channels, primarily because this reduces the extent of branching and 
neighbourhood touching - two major causes of propagation leakage - in much the same 

way as erosion. 

This pre-processing is performed by first thresholding the image, then computing a 
5 distance map using a Chamfer algorithm, and finally thresholding the distance map. 
Using the Chamfer algorithm eliminates the need to administer multiple passes of FMM, 
as proposed by Deschamps et al and the additional complexity that it brings 
(T. Deschamps and L.D. Cohen, "Fast extraction of minimal path in 3D images and 
applications to virtual endoscopy" Medical Image Analysis, 5: 281-299, 2001. 
10 Although the resultant distances are not exactly Euclidean, the approximation was found 
to be adequate as only the order of distances, rather than the distance values per se, are 
of interest to the method of the preferred embodiment, and using any of the common 
Chamfer matrices (e.g. G. Borgefors. Distance transformations in arbitrary dimensions. 
Comput Vision Graphics Image Process, 27: 321-345, 1984) that order is preserved on 
15 the size of grids for CT or MRI scans is suitable. 

Findine a Minimal cost path 

The minimal cost approach is both effective and robust for extracting curvilinear 
structures such as the vasculature and the colon. It is robust to many conditions such as 
vessel stenosis, partial volume effects and noise. Unfortunately, minimizing the path 
20 integral of a cost function is an intractable NP-complete problem. 

Simulated wave propagation has recently emerged as a desirable approach to 
finding the minimal cost path. Not only is it efficient, but it can also give sub-grid 
accuracy. The latter is a distinct advantage over the more traditional graph search 
algorithms as it overcomes an inherent ambiguity associated with a discrete grid. 
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The most commonly used method for efficiently tracking a propagating front is the 
Fast Marching Method (FMM) proposed by Sethian (J.A.Sethian. A fast marching level 
set method for monotcnically advancing fronts. Proc. of the National Academy of 
Sciences of the USA, 93(4):1591-1595). Provided that the speed of the propagation is 

5 monotonic, interface evolution is governed by the Eikonal equation |Vr|F = / , where T 

is the arrival time of the front, and F its speed. The front speed can be defined as: 

G *I 

F= ^ + a, 

max(G a * /) 

where / is the image, G a is a Gaussian kernel and a is a small constant. A small global 

directional bias term 0 can be optionally added to speed up the fast marching process, 
10 but a larger value of j8 impacts the accuracy. One of the purposes of normalising of the 

» 

intensity is to minimise the need to adjust parameters such as a and jS across different 
data sets. 

Efficient entropy-satisfying numerical methods have been developed by Sethian 
and his associates based on the Hyperbolic Conservation Laws. After a forward 

15 marching phase, backtracking from anywhere in the field in the direction of steepest 
decent will reach the origin of the propagation. Thus, from a pre-determined end point, 
a minimal cost path can be found via wave propagation simulated using FMM. 
The safest way to avoid spurious tortuosity being introduced, while ensuring that any 
real tortuosity is not underrepresented by the model, is to segment out any bumps, 

20 branches or foreign objects in the neighbourhood. To be able to recognise those, 
however, one needs to apply an a priori knowledge about the desired and irrelevant 
structures. 
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A tubular shape model: Active-surface-based tubular deformable model can be 
utilised for the recognition. Suitable models can be found in D. Terzopoulos, A. Witkin, 
and M. Kass. Constraints on deformable models: Recovering 3D shape and nonrigid 
motion. Artificial Intelligence, 36(1): 91-123, 1988. The axial constraints of this model 
5 means it will reject many sudden changes in the axial direction. Consequently, the 
irrelevant structures that are filtered out will involve those that include sudden changes 
in diameter. 

Local orientation of the vessel track One commonality among vessels, ducts, 
bronchi and the colon are that they have circular or elliptical cross sections and smoothly 

10 varying radii. To exploit this fact, a tubular shape model is used as a vehicle to carry 
knowledge about the desired structure for the purpose of filtering out any irrelevant 
bumps, branches or foreign objects in the neighbourhood. The knowledge is embedded 
not only in the mesh structure, but also in the internal forces of the model. The 
combination of the model's intra-ring forces (explained below) and the inflation force 

15 alone favours such an outcome that each cross-sectional ring of the resultant tubular 

model has a nearly constant area. Although the inter-ring forces help enforce this, in the 
present application they only play a secondary role. The internal forces can be designed 
so that the constant area favoured by the model approximately corresponds to that of the 
structure being modelled. The use of a deflation force (see below) near the end of the 

20 deformation provides further assurance. The image force localizes the model. 

A feature that makes the model simple and computationally efficient is that the 
model's deformation takes place in a transformed image. As mentioned above, the 
original data are resampled in planes perpendicular to the local orientation of the MCP. 
The resampled data are then stacked up to form a new, transformed data volume. An 

25 initial thin tube is constructed in the middle of the new data and allowed to evolve to 
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minimize the following "energy" functional E of the model in a space of permissible 
deformation: 
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where v(s, r) = (x(s, r), y(s, r), z(s, r)) is a parametric surface on a parameter region 
5 12, s and r are the parameterisation in the cross-sectional-tangential and axial directions 
respectively, P(v) is the potential associated with external forces and can be defined as 

- |V/(v) | , where / is the image. wi 0 , w 0 i, w n , w 2 o, w 0 2 control surface properties of 

tension, rigidity and resistance to twist (they are not necessarily constants). Used in our 
model, the coefficients wio and W20 encode the strengths of the intra-ring forces 
10 mentioned above, while >v 0 i and w 0 2 represent those of the inter-ring forces. w\\ is 
associated with a combination of intra and interring forces. An inflation force is also 
used throughout the process. 

A minimum of E can be reached by solving the associated Euler-Lagrange 
equation 
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where s and r are the parameterisation in the cross-sectional-tangential and axial 
directions respectively. Similar methods can be found in "Finite-Element Methods for 
Active Contour Models and Balloons for 2-D and 3-D Images", by L. Cohen and I. 
Cohen, IEEE Transactions on Pattern Analysis and Machine Intelligence, 15(1 1):1 131- 



20 1147, Nov 1993. 
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Using image gradients as the potential, local minima of the functional are usually 
encountered when the surface passes through some gradients. A technique that can be 
employed to speed up the model's convergence is that, since the value at any position in 
the distance field from a binary image corresponds to the smallest distance from the 
5 point to the background, it is possible to use the distance value at the points on the MCP, 
as the initial radii in the shape model described above. This technique promotes a 
quicker convergence. 

After the model converges at such an energy minimum, the medial axis is 
transformed back to the coordinate system of the original data. This transformed axis is 
10 the curvilinear model that the present method produces. 

This method provides a high degree of accuracy in length, angle and tortuosity 
measurements for certain surgical planning (e.g. endovascular) purposes providing all 
possible measurements concerning a tubular entity, potentially interfering irrelevancies 
can be filtered out before extracting a model suitable for the measurements. The method 
15 provides an implementation of the curvilinear approach and demonstrates that it is 

resistant to introducing spurious curvatures, those that are not due to any genuine change 
in the local orientation of the object of interest, while faithfully reproducing "real" high 
curvatures, those that are actually part of the object in question. 

Further improvements can be made. In the implementation presented above, a 
20 deformable tubular model and a front propagation approach are the two main elements 
used in the filtering and the extraction. However, it is possible to modify or replace 
either of these. For example, incorporating adaptivity in w 0/ in equation 1, in a fashion 

such as that expressed in the equation below can help improve the robustness of this 
approach. 
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vv 0/ (i ( i + i)+a[max(a.,a. + 7)- min(a,,a (+/ )]+/?, (3) 



where a 2 - and a /+/ are areas of two successive rings of the model. However, experience 
has suggested that this adaptivity is not strictly necessary. 

An a Priori estimation of the vessel should be used to give such a scale of smoothing 

5 that it facilitates real edge finding, such as a = — , where d ce is the estimated width of 

2 

the propagating channel. An inflation force adaptive to the sign and magnitude of the 

gradient can also be used. 

The preferred embodiment is preferably implemented on a standard workstation 
type computer arrangement utilising standard computer graphics programming 

10 languages. A suitable hardware arrangement for programming can be as shown in Fig. 4 
with central processing unit 41, memory 42, Disk store 43 and IO Device Controller 44 
arranged around a central bus 45. The arrangement 40 can be suitably programmed in 
the usual manner to carry out the operations of the preferred embdodiment. 

The foregoing describes only preferred embodiments of the present invention. 

15 Modifications, obvious to those skilled in the art, can be made thereto without departing 
from the scope of the invention. 
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THE CLAIMS DEFINING THE INVENTION ARE AS FOLLOWS :- 

1 . A method of constructing a curvilinear model of a tubular structure from a 3- 
dimensional data set, the method comprising the steps of: 

(a) forming a substantially minimal cost path for the tubular structure by the 
5 steps of: 

(al) constructing a distance field; 

(a2) utilising said distance field to construct a propagation map; 
(a3) utilising said propagation map to construct said substantially 
minimal cost path; 

10 (b) for each point on the substantially minimal cost path, providing the steps of: 

(bl) sampling said 3-dimensional data set in a substantially locally 
perpendicular axis to the minimal cost path; 

(b2) constructing a deformable tube around the substantially minimal cost 
path determined by alterations in pixel values around said minimal cost path; 
15 (b3) extracting an axial measure of said deformable tube as said 

curvilinear model of said tubular structure. 

2. A method as claimed in claim 1 further comprising the step (b4) of transforming 
the curvilinear model back into the original data set. 

3.. A method as claimed in claim 1 wherein said step (b2) further comprises the 
20 steps of constructing a resample the data column around the substantially minimal cost 
path and then constructing said deformable tube within said resampled data column . 

4. A method as claimed in claim 1 wherein said distance field is constructed 
utilising a Chamfer algorithm for a distance measure. 
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5. A method as claimed in claim 4 wherein said distance field is thresholded. 

6. A method as claimed in claim 1 wherein said propagation map is constructed 
utilising Fast Marching Method techniques. 

7. A method as claimed in any previous claim 1 wherein said tubular structure 

5 comprises a tubular structure within the human body and said 3 -dimensional data set 
comprises a scan of a portion of the human body. 

8. An apparatus programmed to carry out the operations of any of claims 1 to 7. 
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